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The  purpose  of  this  study  is  to  present  a  2D  transient  numerical  model  to  predict  the  dynamic  behavior  of 
a  tubular  SOFC.  In  this  model,  the  transient  conservation  equations  (momentum,  species  and  energy  equa¬ 
tions)  are  solved  numerically  and  electrical  and  electrochemical  outputs  are  calculated  with  an  equivalent 
electrical  circuit  for  the  cell.  The  developed  model  determines  the  cell  electrical  and  thermal  responses 
to  the  variation  of  load  current.  Also  it  predicts  the  local  EMF,  state  variables  (pressure,  temperature  and 
species  concentration)  and  cell  performance  for  different  cell  load  currents.  Using  this  comprehensive 
model  the  dynamic  behavior  of  Tubular  SOFC  is  studied.  First  an  initial  steady  state  operating  condition 
is  set  for  the  SOFC  model  and  then  the  time  response  of  the  fuel  cell  to  changes  of  some  interested  input 
parameters  (like  electrical  load)  is  analyzed.  The  simulation  starts  when  the  cell  is  at  the  steady  state  in 
a  specific  output  load.  When  the  load  step  change  takes  place,  the  solution  continues  to  reach  to  the  new 
steady  state  condition.  Then  the  cell  transient  behavior  is  analyzed.  The  results  show  that  when  the  load 
current  is  stepped  up,  the  output  voltage  decreases  to  a  new  steady  state  voltage  in  about  67  min. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Solid  oxide  fuel  cells  (SOFCs)  are  a  highly  efficient,  environmen¬ 
tally  benign  method  of  electric  power  production.  SOFC  systems 
have  been  proposed  for  electric  utility  power  generation  in  both 
large  central  station  power  plants  and  distributed  generation 
stations  [1].  Considerable  amounts  of  research  work  have  been 
conducted  to  SOFC,  making  the  SOFC  close  to  the  commercial  appli¬ 
cations.  One  of  the  challenges  for  application  of  SOFCs  is  their 
relatively  slow  response  to  the  input  parameters  time  variation. 
Understanding  the  transient  behavior  of  SOFCs  is  also  important 
for  control  of  stationary  utility  generators  during  power  system 
faults,  surges,  and  switching.  Therefore,  the  transient  modeling  of 
the  fuel  cell  is  useful  to  predict  the  cell  dynamic  behavior  and  oper¬ 
ation.  In  1994,  Achenbach  [2]  analyzed  the  dynamic  operation  of  a 
planar  solid  oxide  fuel  cell.  He  examined  the  transient  cell  voltage 
performance  due  to  temperature  changes  and  current  density  with 
lumped  assumption  for  the  cell  temperature  distribution.  Hall  and 
Colclaser  [3]  also  developed  a  thermodynamic  model  for  prediction 
of  transient  operation  of  the  Tubular  SOFC.  Sedghisigarchi  and  Feli- 
achi  [4]  combined  heat  transfer  dynamics  and  species  dynamics  to 
form  a  new  dynamic  model.  Also  Xue  et  al.  [5],  considered  a  one¬ 
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dimensional  transient  model  for  heat  and  mass  transfer  simulation 
assuming  an  electrical  circuit  include  the  Ohmic  resistances  and 
capacitors  for  the  energy  storage  mode  of  operation.  Iora  et  al.  [6] 
considered  the  internal  reforming/shifting  reactions  in  fuel  channel 
in  their  study.  Qi  et  al.  [7,8]  developed  a  quasi  2D  model  of  a  tubu¬ 
lar  SOFC  based  on  the  changes  along  the  gas  flow  direction  using 
the  control  volume  (CV)  approach.  In  their  model,  the  cell  length 
is  divided  to  several  serial  segments.  Each  segments  includes  five 
CVs,  i.e.  air  tube,  air  channel  (in  two  sections),  cell  body  (electrolyte 
and  electrodes)  and  fuel  channel.  They  obtained  a  non-linear  set 
of  differential  equations  for  the  heat  and  mass  transfer  as  well  as 
electrical  and  electrochemical  variables.  By  solving  these  equations 
simultaneously,  they  calculated  the  cell  time  response  to  the  load 
change.  One  obvious  weakness  of  these  dynamic  models  is  that  they 
used  constant  heat  and  mass  transfer  coefficients  based  on  a  fully 
developed  flow  approximation  at  constant  wall  temperature  and 
mass  flux.  In  fact  the  tubular  SOFC  is  a  heat-generating  tube  with 
different  flow  streams  on  both  the  inner  and  outer  side.  We  will 
avoid  the  fully  developed  flow  approximation  of  Nusselt  and  Sher¬ 
wood  numbers  at  constant  wall  temperature/concentration  in  this 
2D  transient  model  by  a  complete  field  solution  of  the  governing 
equations  for  the  heat  and  mass  transfer  in  the  entire  domain  of  a 
tubular  SOFC  working  in  a  cell  stack.  Also  the  radial  gradient  of  the 
state  variables  (i.e.  temperature  and  species  mole  fractions)  in  air 
and  fuel  channels  is  predicted  and  therefore  the  time  they  need  to 
diffuse  in  whole  domain  is  taken  in  to  account  in  present  dynamic 
model. 
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Nomenclature 

A  surface  area  (m2) 

Cp  specific  heat  capacity  (J  kg-1  k-1 ) 

Dj  m  diffusion  coefficient  of  jth  species  into  the  left  gases 

of  a  mixture  (m2  s-1 ) 

E  cell  ideal  voltage  (V) 

Ea  activation  energy  of  anode  (J  mol-1 ) 

Ec  activation  energy  of  cathode  (J  mol-1 ) 

F  Faraday’s  constant  (C  mol-1 ) 

G  Gibbs  function  (J  kg-1 ) 

G°  standard  Gibbs  function  (J  kg-1 ) 

T  temperature  (I<) 

I  current  (A) 

i  current  density  (A  m-2 ) 

k  pre-exponential  factor  (A  nrr 2 ) 

fn  mass  flux  (kg  m-2  s-1) 

M  molecular  mass  (g  mol-1) 

ne  number  of  electrons  participating  in  electrochemi¬ 

cal  reaction 
P  pressure  (Pa) 

Q.  volumetric  heat  source  ( W  nrr3 ) 

q  heat  source  (W) 

r  radius  (m) 

R  electrical  resistance  (£2) 

U  utilization  percentage  (0-1 ) 

u  axial  velocity  (m  s-1 ) 

v  radial  velocity  (m  s-1 ) 

V  voltage  (V) 

Xj  jth  species  mol  fraction 

Yj  jth  species  mass  fraction 

Greek  symbols 

I I  electrical  losses  (V) 

8  thickness  of  layers  (anode,  cathode  and  electrolyte) 

(m) 

pp  electrodes  and  electrolyte  electrical  resistivity  in 
nodep  (£2m) 

0  cell  peripheral  angel  (rad) 

p  density  (kg nr3) 

jjl  viscosity  (Pas) 

A  thermal  conductivity  (W  nr1  K-1 ) 

AH  hydrogen  heat  value  (j  mol-1 ) 

91  universal  gas  constant  (J  mol-1  K-1 ) 

A  V  volume  of  an  element  in  the  electrodes  or  electrolyte 

(m3) 

Subscripts 
a  anode 

c  cathode 

e  electrolyte 

act  related  to  activation  loss 

p  node  p  in  the  network  circuit 

£2  related  to  Ohmic  loss 


2.  Model  description 

2.1  Geometrical  model 

The  computational  domain  and  the  flow  streams  of  fuel  and 
oxidant  in  a  typical  tubular  SOFC  is  shown  in  Fig.  1.  Due  to  symme¬ 
try,  only  half  of  the  cell  unit  (between  the  cell  symmetry  axis  and 
the  symmetry  line  between  two  adjacent  cells)  is  being  analyzed. 
Fig.  2  shows  the  detailed  view  of  the  axi-symmetric  computational 


Fig.  1.  The  schematic  diagram  of  a  tubular  SOFC  and  the  computational  domain. 


domain  and  the  cell  overall  heat  and  mass  flows.  As  shown  by  the 
Figure,  the  computational  domain  includes  the  air  and  fuel  chan¬ 
nels,  anode,  cathode,  electrolyte  and  the  supporting  tube  layers. 
The  fuel  entrance  is  at  the  closed  end  of  the  cell.  Table  1  presents 
the  geometrical  characteristics  of  the  tubular  SOFC  [9,10].  The  gov¬ 
erning  electrochemical,  electrical  as  well  as  heat  and  mass  transfer 
equations  are  as  follows: 

2.2.  Electrochemical  model 


The  local  ideal  electromotive  force  (EMF)  for  the  hydrogen  elec¬ 
trochemical  reaction  is  given  by  [12]: 

-AG  -AC0  9tT  (Po2/P0)°'5(Ph2/P0) 

^  =  ~rieFr  =  2F  ^  2F  n - P^Jpo -  (1) 

where  P°  is  the  atmospheric  pressure.  Activation  losses  in  anode 
and  cathode  of  an  elemental  area  ( AA)  are  as  follows,  respectively 
[2]: 


he  = 

ha  = 


I 

A/\c((4F/fiP)fcc(Po2/P)0-25  exp {-EclRT)) 

I 

A/\a((4P/Rr)ka(PH2/P)0-25exp(-£a/Rr)) 


=  Kctl  (2) 

=  Ractl  (3) 


where  P  and  T  are  local  pressure  and  temperature,  respectively.  The 
hydrogen  partial  pressure  is  used  in  fuel  channel  while  the  oxygen 
partial  pressure  is  considered  in  the  air  channel.  The  coefficients  of 
Eqs.  (2)  and  (3)  are  given  in  Table  2.  The  concentration  polarization 
due  to  the  non-uniform  distribution  of  the  species  existing  in  the 
air  and  fuel  flow  passages  is  automatically  taken  into  considera¬ 
tion  by  the  mass  transfer  model,  but  the  concentration  polarization 
related  to  the  diffusion  of  the  species  within  the  porous  electrodes 
is  ignored  [3]. 


Table  1 

Geometrical  characteristics  of  the  tubular  SOFC  unit  [9]. 


Thickness  (|jim) 

Supporting  tube 

1500 

Cathode 

1000 

Electrolyte 

50 

Anode 

150 

Diameter  (mm) 

Inner  side  of  air-inducing  tube 

8.0 

Outer  side  of  air-inducing  tube 

9.0 

Inner  side  of  supporting  tube 

13.8 

Outer  side  of  Anode 

19.2 

Outer  diameter  of  fuel  channel 

29.2 

Length  (mm) 

Tubular  cell  unit  500 
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Fig.  2.  Schematic  of  flow  field  in  the  computational  domain. 


Table  2 

Coefficients  of  the  activation  polarization  equation  [11  ]. 


Coefficients 

values 

kc 

1.49  x  1010  Am-2 

ka 

2.13  x  108  Am-2 

Ec 

1.6  x  105  J  mol-1 

Ea 

1.1  x  105  Jmol-1 

noted  that  the  interconnect  resistances  was  neglected  in  compare 
with  other  resistances  in  the  single  cell.  In  order  to  calculate  the  cur¬ 
rent  and  voltage  distribution  in  a  specified  load  current,  Kirchhoff  s 
law  is  applied  to  each  node.  For  each  node,  an  equation  associating 
the  potential  of  the  central  grid  point  P  in  anode  with  the  poten¬ 
tials  of  its  neighboring  points  (east,  west,  north,  and  south)  and 
also  the  corresponding  grid  P  in  the  cathode  is  obtained.  With  all 
of  the  equations  for  the  descritized  grids  in  both  the  cathode  and 
anode  given,  a  matrix  equation  for  potentials  is  created.  The  voltage 
in  cathode  connector  is  set  to  zero  (as  a  reference  voltage)  and  the 
total  current  taken  out  from  the  cell  is  prescribed.  Then  the  poten¬ 
tials  in  each  node  and  also  in  anode  connector  are  calculated  by 
solving  the  matrix  of  potential  equations. 

The  electrical  resistances  in  radial  direction  for  the  anode,  cath¬ 
ode  and  electrolyte  are,  respectively: 

” p  —  *act 

( _ 1 _ 

V  AA,((4F/RTMPh2/P)0-25  exp(-EalRT)) 


Fig.  3.  Network  circuit  for  a  discretized  element  of  the  tubular  SOFC. 


2.3.  Electrical  model 

Fig.  3  shows  the  equivalent  network  circuit  for  a  discretized  ele¬ 
ment  of  the  tubular  SOFC.  The  equivalent  electrical  circuit  is  used  to 
calculate  the  cell  current  and  voltage  distribution  in  a  specified  load 
current.  The  circuit  includes  a  set  of  electrical  resistances  and  local 
electrical  power  sources.  The  schematic  view  of  the  overall  equiv¬ 
alent  circuit  for  half  of  the  cell  is  depicted  in  Fig.  4.  It  should  be 


e 


Fig.  4.  Schematic  view  of  the  overall  equivalent  circuit  for  half  of  the  cell. 


RCP  =  Rh 


+  R 


C 

act 


(  A/4c((4F/Rr)kc(Po2/P)0'25 


exp  (-EJRT)) 


(5) 


(6) 


A/4  is  the  average  peripheral  area  of  each  element  of  electrodes  or 
electrolyte  and  is  AA  =  (rin  +  rout )  x  d Ox  dx/2  (Fig.  3);  where  rin  and 
Tout  are  inner  and  outer  radii  of  each  layer,  respectively.  pv  is  elec¬ 
trical  resistivity  of  the  cell  material  in  node  p,  which  is  constant  for 
the  cathode  and  anode,  but  is  temperature-dependent  for  the  elec¬ 
trolyte,  as  it  is  given  in  Table  3  [ll].InEqs.(4)and(5),PH2  andPo2  are 
partial  pressures  of  H2  and  02  species  at  electrode-gas  interfaces, 
respectively  and  T  is  the  corresponding  electrode  temperature. 

The  longitudinal  and  peripheral  resistances  for  the  anode  are 
calculated  by,  respectively: 


RaPE  = 


PpVadO 
8  a  dx 


nci  _ 

nPN  ~ 


_Pp dx_ 

Saradff 


(7) 

(8) 


The  same  set  of  relations  (Eqs.  (7)  and  (8))  is  used  for  the  cathode. 
The  voltage  difference  between  anode  and  cathode  in  any  node  is 
given  by 


Vpa  -Vcp=Ev-  {R“  +Rep  +Rcp)Ip 


(9) 


2.4.  Heat  and  mass  transfer  model 


To  couple  the  flow  and  heat/mass  transfer  fields,  a  two- 
dimensional  axi-symmetric  computational  domain  is  created  (as 


Table  3 

Solid  components  properties  of  the  cell  [11]. 
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Cell  component 

Thermal  conductivity  (W(m  K)-1 ) 

Heat  capacity  (J  (kg K)-1 ) 

Density  (kg  m-3) 

Electrical  resistivity  (£2  m) 

Cathode 

6.0 

623.0 

4930.0 

0.0186 

Electrolyte 

2.7 

623.0 

5710.0 

0.3685  +  0.002838  exp(  10300 /T) 

Anode 

11.0 

623.0 

4460.0 

0.0014 

Supporting  tube 

1.1 

623.0 

4460.0 

- 

Air  inducing  tube 

1.1 

623.0 

4320.0 

- 

shown  in  Fig.  2).  Laminar  flow  assumption  is  used  for  the  heat  and 
mass  transfer  in  air  and  fuel  channels.  This  assumption  is  due  to 
rather  low  flow  velocity  and  Reynolds  number  (maximum  1000 
in  air  channel)  in  both  air  and  fuel  channels  [10,15].  Conservation 
equations  are  applied  universally  to  entire  computational  domain, 
however  zero  velocities  is  assigned  to  solid  area  in  the  numeri¬ 
cal  treatment.  The  transient  momentum,  heat  and  mass  transfer 
governing  equations  in  cylindrical  coordinate  are  used  as  follows, 
respectively: 


In  anode-electrolyte  interface:  the  thermodynamic  heat  source 
is  due  to  hydrogen  electrochemical  reaction  [12]: 

flTh=(^-£n,p)/P  (17) 

En,p  is  the  local  Nernst  voltage  and  Ip  is  the  local  current  in  node 
P- 

The  mass  flux  into/out  of  the  fuel/anode  boundary  as  well  as 
air/cathode  boundary  is  as  follows,  respectively: 


d{pu)  d{puu)  1  d(rpvu) 
dt  +  dx  +  r  dr 


mu2  ~  ™h2o  =  Mh2  —  -  Mh2o  2f 


(18) 


dp 


du 


dx  +  dx  dx  )  +  r  dr  \  ^  dr 


1  d 


du 


du 


1  d 


+  3x  dx  j  +  r  3r  l  m  dx 


dv 


(10) 


d{pv)  (  d{puv )  (  1  d(rpvv) 
~d F  +  ~dx~  +  7  dr 


dp 


dv 


dr  +  dx[fldx  +rdr[mdr 


1  3 


dv 


du 


1  d 


+  37  ^37  +  ~rdr  \^Tr  -  r 


dv 


2  /IV 


(11) 


d(pCpT)  d(pCpuT)  1  d(rpCpVT) 

dt  dx  r  dr 


d  (  dT\  1  d  dT 
dx  V  dx)  +  r  dr  V  r  dr 


+  Q 


(12) 


d(puYj)  1  d(rpvYj) 
dt  +  dx  +  r  dr 

d  /  dYj\  U  /  dYj 
~  (  PDfm  ^  )  +  r  (  rPDj,m 


9x 


r  9r 


(13) 


In  Eq.  (12)  the  source  term  (Q.)  is  treated  as  follows: 

In  anode  and  cathode:  the  Joule  volumetric  heating  source  term 
due  to  Ohmic  resistance  and  activation  polarization  is  [11  ]: 


f)0  P  P 

-  AVH 

(14) 

RZ  ■  \l 

II 

< 

3n  ^ 

(15) 

In  electrolyte:  the  Joule  volumetric  heating  source  term  is  due 
to  Ohmic  resistance: 


Rl_ 1 

AVI 


(16) 


m  o2=Mo2^p  (19) 

2.4.1.  Physical  and  other  properties 

All  transport  and  thermo-physical  properties  p,  /x,  Cp ,  A  and  Dy  of 
the  fuel  and  air  are  calculated  as  functions  of  the  local  temperature, 
pressure,  and  species  concentrations.  All  the  gaseous  components 
are  assumed  to  behave  as  perfect  gases,  and  the  local  properties  of 
the  gas  mixture  are  calculated  according  to  the  relations  given  in 
reference  [16].  Table  3  gives  the  values  used  for  heat  conductivity, 
density  and  heat  capacity  and  also  the  ionic  and  electric  resistivities 
for  the  solid  components  [11  ]. 


2.4.2.  Boundary  conditions 

The  general  boundary  conditions  for  the  computational  domain 
(Fig.  2)  are: 

The  inlet  velocity  of  the  fuel  and  air  flow  for  a  specified  utility 
factor  is  given  by  [11]: 


_  f  A cellicell  A  NTfuel 

fuel~  \2FUH2XH2Afuel )  Pfuel 


(20) 


^air 


AcelPcell  |  ^Tair 

4FUo2Xo2Aair  J  Pa[r 


(21) 


The  species  mass  flux  in  air/cathode  and  fuel/anode  interfaces 
(due  to  the  local  convection  velocity  and  species  mass  diffusion)  is 
as  follows,  respectively  [11]: 

mf  -  Dhairpf  ^  +  pairY}vair  (22) 

rfue,  =  -DhmrfeldJ+^YjVf  (23) 

where 


and 


Vfuel  ~ 


E  rf* 


(24) 


(25) 


where  AVP  is  the  volume  of  an  element  in  the  electrodes  and  elec-  The  fuel  exterior  boundary  (the  boundary  between  adjacent 

trolyte  and  is  A  Vp  =  A Ap  x  8  and  8  is  the  thickness  of  each  layers.  cells)  is  assumed  to  be  symmetrical. 
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Fig.  5.  Flowchart  of  the  cell  transient  simulation  process. 


3.  Simulation  procedure  of  the  cell  transient  operation 

The  flowchart  for  the  overall  simulation  process  of  the  cell  tran¬ 
sient  operation  is  shown  in  Fig.  5.  The  simulation  is  based  on 
a  finite  volume  based  code  that  uses  separate  modules  for  each 
sub-model  (electrical,  electrochemical  and  heat  and  mass  transfer 
sub-models).  The  steady  state  solution  of  the  cell  is  set  as  the  initial 
value  for  the  transient  analysis.  As  shown  in  Fig.  5,  at  each  time  step 
an  internal  iteration  is  performed  and  continues  until  it  converge. 
Then  the  code  proceeds  to  the  next  time  step  utilizing  previous 
results. 

The  governing  transient  transport  equations  (momentum, 
energy  and  species  conservation  equations)  are  solved  by  the  SIM¬ 
PLE  algorithm.  In  this  method  the  time  steps  do  not  affect  the 
stability  of  the  solution;  however,  it  affects  the  accuracy  of  the 
results  [  13  ].  The  computational  domain  is  discretized  by  a  220  x  640 
grid  in  radial  and  longitudinal  directions,  respectively.  The  grid 
independency  of  the  results  is  examined  and  insured.  The  grids 
are  uniform  in  longitudinal  direction  and  non-uniform  in  the  radial 
direction.  This  is  due  to  very  low  thickness  of  the  anode  and  elec¬ 
trolyte  layers  relative  to  the  air  and  fuel  channel  width. 


4.  Results  and  discussion 


the  experimental  data  especially  for  current  density  greater  than 
250  mA  cm-2.  Results  for  current  density  less  than  250  mA cm-2 
show  about  12%  deviation  with  experimental  data.  This  is  due  to 
the  assumption  and  formulation  of  the  activation  and  concentration 
polarization  as  well  as  the  accuracy  of  numerical  solution  proce¬ 
dure.  Furthermore,  the  present  model  results  show  less  deviation 
from  the  experimental  results  than  Li  et  al.  results  (about  5%  more 
accurate).  This  is  due  to  different  treatment  of  the  activation  and 
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For  verification  purpose,  the  current  study  polarization  curve 
is  compared  to  experimental  data  of  Flagiwara  et  al.  [14]  and  the 
steady  state  model  of  Li  et  al.  [9]  in  Fig.  6.  As  shown  the  obtained 
results  show  good  agreement,  less  than  0.25%  deviation,  with 
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Fig.  6.  Cell  voltage  versus  current  density. 
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Table  4 

The  cell  operating  conditions  in  steady  state  model  [9]. 
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Operating  parameters 

Values 

The  cell  current  density/output  current 

3500Am_2/104  A 

Input  air  temperature 

600  °C 

Input  fuel  temperature 

840  °C 

Air  and  fuel  pressure 

1.013  atm 

Input  air  composition 

98%  H2,  2%  H20 

Hydrogen  utility  factor  in  fuel  channel/fuel  mass  flow  rate 

0.85/0.00044  gs-1 

Oxygen  utility  factor  in  air  channel/air  mass  flow  rate 

0.167/0.054  gs-1 

concentration  losses  in  two  models.  They  assumed  the  sum  of  the 
activation  polarization  and  the  concentration  polarization  occur¬ 
ring  inside  the  porous  electrodes  is  constant  and  is  proportional 
to  the  electric  current  density.  In  general,  it  is  well  documented 
that  sophisticated  formulation  and  treatment  of  the  activation  and 
concentration  loss  causes  more  accurate  simulation  results. 

4.1.  Steady  state  results 

The  steady  state  operation  results  of  the  developed  model  are 
presented  in  this  section.  Input  air  and  fuel  conditions  and  also 
other  setups  of  the  cell  for  the  steady  state  model  are  presented 
in  Table  4.  These  steady  state  results  are  used  as  the  initial  values 
for  the  transient  simulation. 

Fig.  7  shows  the  cell  temperature  contours  for  load  current 
density  of  3500  Am-2.  The  flow  temperature  rises  downstream  of 
the  inlet  and  approaches  to  about  1000  °C  due  to  the  heat  that  is 
generated  in  the  solid  part  and  is  conducted  into  the  flow.  Then 
temperature  slightly  decreases  toward  the  cell  closed  end.  Because 
the  air  that  absorbs  the  major  portion  of  the  generated  heat  because 
of  its  higher  mass  flow  rate.  Therefore  the  radial  temperature  gra¬ 
dient  in  the  air  flow  is  more  then  the  fuel  flow.  In  the  closed  end 
region  of  the  cell,  non-uniformity  of  air  temperature  is  reduced  due 
to  the  mixing  related  to  the  flow  recirculation  there. 

The  relative  pressure  contours  in  air  channel  is  depicted  in  Fig.  8. 
As  shown  by  figure  the  pressure  loss  in  air  channel  is  calculated  to 
be  about  637  Pa.  The  input  air  gage  pressure  is  set  to  1280  Pa. 

Fig.  9  shows  the  oxygen  mass  fraction  contours  in  air  channel. 
Diffusion  of  oxygen  from  the  core  region  to  the  support  tube  surface 
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Fig.  10.  Hydrogen  and  vapor  mass  fraction  contours  in  the  fuel  channel. 

forms  a  significant  radial  mass  fraction  gradient.  As  shown  by  the 
figure  at  the  air  tube  outlet  the  oxygen  mass  fraction  reaches  to 
0.207. 

Fig.  10  shows  the  hydrogen  and  vapor  mass  fraction  contours  in 
the  fuel  channel.  In  the  axial  direction,  the  mass  fraction  of  water 
vapor  increases  along  the  fuel  channel  from  0.13  up  to  0.98,  which 
is  due  to  strong  production  rate  of  vapor  from  the  anode  wall. 
In  the  radial  direction,  water  vapor  and  hydrogen  concentrations 
do  not  show  significant  gradients.  This  is  due  to  the  strong  diffu¬ 
sion  of  hydrogen  in  the  fuel  stream  and  very  low  fuel  mass  flow 
rate. 

The  flow  streamlines  in  fuel  and  air  channels  are  shown  in  Fig.  11 . 
As  shown  the  injected  mass  from  the  anode  wall  is  comparable  with 
the  main  mass  flow  in  fuel  channel.  This  is  because  of  low  mass  flow 
rate  of  input  fuel  (due  to  high  utility  factor  of  hydrogen).  However, 
in  the  air  channel  the  consumed  mass  by  the  cathode  wall  is  very 
low  compared  with  the  input  air  mass  flow  rate. 

Fig.  12  shows  the  axial  distribution  of  the  cell  local  Nernst  volt¬ 
age  and  current  density  in  the  average  load  current  density  of 
3500  Am-2.  Maximum  Nernst  voltage  is  about  1.13  V  at  the  closed 
end  of  the  cell  where  the  fuel  is  coming  in.  in  the  inlet  of  the  fuel 
the  H2  mass  fraction  is  maximum  and  therefore  Nernst  potential  is 
maximum  there.  Also  current  density  decreases  along  the  cell  by 
decreasing  the  cell  local  voltage. 

Generally  the  steady  state  results  obtained  from  present 
dynamic  model  are  in  good  agreement  with  the  steady  state  models 
previously  reported  by  Li  et  al.  and  also  experimental  data  published 
by  Hagiwara  et  al.  [9,14]. 


0  10  20  30  40  50 

x  (cm) 

Fig.  7.  Cell  temperature  contours  (°C)  in  the  whole  computational  domain. 


1120.01  1077.51 

C 

898.876  898.876 

792.047 

634.9^2 

1110.19  1124.35 

1155.9 

1186.48 

1216.67 

127 

1.03 

0  10 

20 

30 

40 

50 

x  (cm) 

Fig.  8.  Relative  pressure  contours  (Pa)  in  the  air  channel  with  the  input  air  gage 
pressure  of  1280  Pa. 
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4.2.  Load  change  results 

The  simulation  was  conducted  to  determine  the  transient 
response  of  the  cell  to  electrical  load  changes.  The  cell  was  subjected 
to  a  step  increase  in  current  density  (from  1500  to  3500 Am-2) 
with  the  fuel  and  air  utility  factors  held  constant.  This  means  that 
the  input  fuel  and  air  mass  flow  rates  increases  with  load  current 
increase. 

Fig.  13  shows  the  axial  temperature  profiles  of  the  cell  during 
the  transient  operation  up  to  steady  state,  for  a  step  electrical  load 
increase  (from  1500  to  3500  Am-2).  As  shown  in  figure,  the  tem¬ 
perature  decreases  at  closed  end  side  of  the  cell  (x  <  0.18  m)  with  an 
increase  of  load  current  but  it  increases  at  the  air  inlet  side;  how¬ 
ever  the  maximum  temperature  decreases  to  about  1000  °C.  Also 
the  cell  maximum  temperature  position  shifts  to  the  right  hand  side 


Fig.  9.  Oxygen  mass  fraction  contours  in  the  air  channel. 


Fig.  11.  The  flow  stream  lines  in  fuel  and  air  channels. 
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Fig.  12.  Axial  distribution  of  local  Nernst  voltage  and  current  density  in  the  load 
current  density  of  3500  A  m- 2 . 


from  0.1  to  0.3  m  (toward  the  air  inlet  throat).  The  reason  is  that  the 
increase  of  the  cell  load  current  leads  to  increase  of  the  air  mass 
flow  rate  (about  130%)  and  its  cooling  effect,  despite  the  increase 
of  the  electrochemical  and  electrical  heat  sources.  Also  maximum 
air-cooling  effect  takes  place  at  the  closed  end  of  the  cell  and  it 
absorbs  maximum  heat  value  from  this  region.  However  at  the  air 
inlet  side,  the  increase  of  heat  sources  causes  the  increase  of  the 
cell  body  temperature  because  the  air-cooling  effect  is  less  than 
thermal  and  electrical  heating  in  this  region. 

Fig.  14  shows  the  cell  closed  end  temperature  history,  under  step 
external  load  change.  As  it  can  be  seen  in  the  figure  the  temperature 
decreases  about  110  °C  by  increase  of  the  load  current.  Also  the  time 
the  cell  temperature  needs  to  reach  to  the  new  steady  state  condi¬ 
tion  (in  current  density  of  3500  Am-2)  is  calculated  about  7000s. 
However,  the  temperature  response  time  obtained  in  present  study 
is  more  than  Hall  et  al.  work  (4500s)  [3].  It  is  due  to  the  radial 
gradient  of  the  state  variables  (i.e.  temperature  and  species  mole 
fraction)  and  the  time  they  need  to  diffuse  in  whole  domain.  Also  an 
over  prediction  of  the  solid-gas  heat  transfer  coefficient  can  induce 
some  errors  in  lumped  or  one-dimensional  models. 

The  electrical  transient  response  of  the  cell  to  the  change  in  load 
current  is  depicted  in  Figs.  15-17.  Fig.  15  represents  the  time  history 
of  the  cell  Nernst  potential.  As  the  external  load  current  goes  up 
from  1500  to  3500  Am-2,  the  cell  Nernst  potential  goes  down  from 
0.934  to  0.926  V.  it  shows  a  fast  response  part  (about  120  s)  initially 
and  then  go  back  to  new  steady  state  value  of  0.926  V  during  a  slow 
transient  response  part  (about  1  h  in  the  current  simulation).  The 
transient  response  in  the  Nernst  potential  is  very  complex,  which 


Fig.  14.  The  cell  closed  end  temperature  history  under  step  external  load  change. 


Fig.  15.  Time  history  of  the  cell  Nernst  potential. 


is  related  to  all  parametric  changes  associated  with  the  entire  fuel 
cell  activity. 

Fig.  16  shows  the  cell  output  voltage  with  time.  After  the  load 
current  experiences  sudden  change,  the  external  load  voltage  goes 
down  from  0.84  to  0.7  V  initially  and  shows  a  very  fast  initial 
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Fig.  13.  The  axial  temperature  profiles  of  the  cell  during  transient  operation. 


Fig.  16.  The  cell  output  voltage  transient  response  to  the  change  in  load  current. 
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Fig.  17.  The  efficiency  and  power  density  history  of  cell  under  step  external  load 
change. 

response.  Then  the  voltage  decays  to  the  new  steady  state  value 
(0.67  V)  gradually,  after  about  67  min.  this  means  that  about  82% 
of  voltage  loss  takes  place  in  a  very  short  time  and  other  remain 
18%,  decreases  during  a  very  long  time,  about  67  min.  it  should 
be  noted  that  because  of  very  slow  thermal  response  of  the  cell, 
the  cell  temperature  distribution  remains  the  same  as  before,  at 
the  first  moment  of  the  load  change.  Therefore  the  output  voltage 
changes  with  the  last  temperature  and  species  conditions  initially. 
Then  the  temperature  and  species  distribution  reaches  to  the  new 
steady  condition  associated  with  the  entire  fuel  cell  heat  and  mass 
sources  and  they  affect  the  electrical  parameters  during  a  long  time 
response  (Fig.  16). 

Fig.  17  depicts  the  efficiency  and  power  density  history  of  the 
fuel  cell  under  step  external  load  change.  As  shown  in  figure,  the 
cell  power  density  increases  from  1250  to  2344  Wrrr2  when  the 
external  load  current  goes  to  3500  Am-2.  It  shows  an  overshoot  of 
about  110%  initially  and  then  it  decays  back  to  the  new  steady  state 
value.  Also  the  cell  efficiency  decreases  from  65%  to  about  52%  after 
load  change. 

5.  Conclusions 

The  newly  developed  2D  transient  model  provides  the  SOFC 
dynamic  response  to  the  cell  electrical  load  change  as  well  as  the 
steady  state  operation.  The  conclusions  drawn  are  as  follows: 

•  The  presented  2D  model  gives  more  realistic  results  than  the 
previous  one-dimensional  and  lumped  models  for  the  dynamic 
operation  of  the  tubular  SOFC,  due  to  considering  the  r-direction 
gradients  of  the  state  variables  (i.e.  temperature  and  species  mole 
fraction). 


•  When  the  load  current  is  stepped  up  from  1500  to  3500  Am-2, 
the  output  voltage  decreases  to  0.67  V  after  about  67  min. 

•  In  a  constant  fuel  and  air  utility  factors,  the  cell  maximum  temper¬ 
ature  decreases  with  increase  of  load  current  and  it  shifts  toward 
the  air  inlet  side. 

•  The  inlet  air  mass  flow  rate  and  thermal  and  electrical  heat 
sources  are  main  parameters  that  determine  the  cell  tempera¬ 
ture  distribution.  The  air-cooling  effect  in  closed  end  of  the  cell  is 
more  than  other  places. 

•  The  temperature  and  species  concentration  time  response  of  the 
cell  is  very  slow.  It  is  calculated  about  7000  s  in  studied  case. 

•  Generally  82%  of  the  cell  electrical  parameters  variations  take 
place  just  after  the  load  change  (very  fast)  and  other  remain  18% 
of  variations  is  during  a  very  long  time  relatively.  This  means  that 
the  effect  of  the  cell  temperature  and  species  concentration  on 
the  electrical  parameters  time  response  is  about  18%  of  the  overall 
their  variations. 

•  In  longitudinal  direction  the  mass  fraction  of  vapor  increases  as 
well  as  the  total  fuel  mass  flow  rate;  while  the  oxygen  mass  frac¬ 
tion  decreases. 

•  Radial  non-uniformity  of  oxygen  concentration  in  the  air  chan¬ 
nel  is  noticeable  so  that  the  related  concentration  polarization  is 
judged  to  be  major.  In  contrast  to  this,  the  radial  non-uniformity 
of  the  hydrogen  concentration  is  not  significant.  Therefore,  con¬ 
centration  polarization  induced  by  the  non-uniform  distribution 
of  species  concentration  in  the  flow  space  is  not  negligible  at  the 
cathode  side  but  it  can  be  ignored  on  the  anode  side. 

•  Pressure  drop  is  about  640  Pa  in  the  air  channel  while  it  is  very 
low  in  fuel  channel  and  can  be  neglected. 
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